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Abstract:  Background:  Methods  for  describing  physics  of  impact  and  ballistics  have  been 
developed  over  a  number  of  decades.  These  include  analytical  mathematical 
representations  as  well  as  modem  computer  simulations. 

Objective:  Recent  and  historic  developments  towards  modeling  of  impact  phenomena 
pertinent  to  tenninal  ballistic  events  are  summarized  and  compared.  Two  classes  of 
physical  problem  are  of  focus:  impact  and  penetration  of  metallic  and/or  ceramic  targets  by 
projectiles,  and  propagation  of  planar  shock  waves  through  solid  material  specimens 
induced  by  collision  with  flyer  plates  or  by  explosive  loading. 


Method:  The  projectile-target  problem  is  analyzed  from  perspectives  of  classical  hydrodynamics,  extensions 
accounting  for  strength,  and  fully  resolved  explicit  dynamics  simulations.  The  planar  impact  test  is  studied  from 
perspectives  of  analytical  solutions  to  Rankine-Hugoniot  equations,  steady  wave  analysis,  and  dynamic  finite 
element  simulations  of  shock  waves  in  material  microstructures.  Key  features  of  each  approach  are  critically 
compared. 


Results:  The  two  classes  of  physical  problem  are  inherently  related  since  material  properties  obtained  from 
analysis  of  the  latter  experiments  are  typical  input  for  models  of  the  former  problem  involving  ballistic 
penetration.  Patents  to  computer  methods  and  ballistic  protection  systems  are  noted. 


Conclusion:  Reduced  order  models  are  shown  to  provide  efficient,  but  often  approximate,  solutions  giving 
insight  into  general  trends.  Modem,  fiilly  resolved  calculations  appear  to  be  the  only  viable  route  to  design  and 
optimization  of  novel  materials  or  stmctures  with  heterogeneous  properties  or  complex  geometries. 


Keywords:  Ballistics,  continuum  mechanics,  shock  physics,  projectiles,  impact,  metals,  ceramics,  crystals. 


1.  INTRODUCTION 

Impact  and  penetration  mechanics  are  physical  phenomena 
that  are  of  high  importance  to  defense  and  industrial 
applications,  including  those  in  mining,  automotive,  and 
aircraft  industries  [1,  2].  Furthermore,  such  phenomena, 
especially  those  involving  hypervelocity  impact,  occur  in 
problems  relevant  to  geology  and  astronomy,  e.g.,  collisions  of 
planetary  objects  or  space  debris  [3].  The  present  work  is 
focused  on  problems  in  terminal  ballistics,  encompassing 
dynamic  interactions  of  projectiles  with  their  intended  targets. 
The  purpose  of  this  paper  is  to  categorize  and  compare  various 
modeling  strategies  that  describe  a  few  characteristic  physical 
problems  in  terminal  ballistics,  noting  the  important  features. 
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relative  advantages  and  disadvantages,  and  key  historic  and 
modem  references  for  each  strategy.  The  intent  is  to  provide  a 
useful  introductory  reference  for  researchers  beginning  work  in 
the  area,  as  well  as  a  source  of  important  equations  and 
references  for  scientists  and  engineers  more  experienced  in  the 
fields  of  ballistics  and  shock  physics.  However,  no  attempt  is 
made  to  cite  all  relevant  works  dealing  with  any  given  topical 
area. 

Relevance  of  this  work  to  recent  patents  in  engineering  is 
as  follows.  This  paper  covers  modeling  techniques— 
formulation  of  governing  differential  equations  and  their 
analytical  or  numerical  solutions— rather  than  development 
of  new  experimental/diagnostic  devices,  munitions,  or 
protection  systems.  Although  the  content  does  not  directly 
invoke  technologies  available  in  recent  patents,  it  does 
supplement  any  such  developments  of  new  technologies.  A 
recent  patent  for  finite  element  computer  methods  used  to 
obtain  numerical  solutions  for  problems  involving  large 
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deformations  as  occurring  under  impact  loading  is  reported 
in  [4],  Recent  patents  for  ballistic  technologies  involving 
armor  ceramics  and  shock  mitigation  include  [5,  6],  Explicit 
listings  of  numerous  other  relevant  patents  for  dynamic 
protection  systems  can  be  found  in  associated  references  to 
industrial  and  aerospace  applications,  for  example  [2,  3], 

The  first  physical  problem  addressed  is  impact  and 
penetration  of  a  solid  target,  typically  metallic  or  ceramic,  by 
a  projectile  such  as  a  bullet  or  rod  at  null  obliquity,  i.e.,  the 
velocity  vector  of  the  the  projectile  is  normal  to  the  plane  of 
the  impacted  surface.  A  comprehensive  book  that  deals  with 
this  problem  and  related  others  is  [7],  including 
experimental,  analytical,  and  numerical  methods  of  ballistics 
research.  Similar  lengthy  references  dealing  with  pertinent 
aspects  include  [8,  9].  In  contrast,  the  present  paper  covers 
only  models— analytical  and  numerical  simulations— with 
experiments  noted  only  in  the  context  of  description  of  the 
physical  problems  and  validation  of  the  analytical  or 
numerical  solutions.  Thus,  this  paper  is  more  concise  and 
focused  than  these  prior  works,  and  it  also  includes  coverage 
of  more  recent  developments,  particularly  in  areas  of 
dimensional  analysis  [10,  11,  12,  13]  and  modem  hydrocode 
simulations  [14,  15,  16].  The  present  description  of 
constitutive  models  of  (poly)crystalline  solids,  of  which 
projectile  and  target  are  comprised,  is  also  more 
sophisticated  than  those  given  in  the  above-mentioned  prior 
works,  accounting  for  geometric  and  material  nonlinearity 
[17]. 

The  second  physical  problem  addressed  is  the  response 
of  a  (poly)crystalline  solid  body  to  dynamic  loading  by  a 
planar  shock  wave.  A  standard  method  for  studying  the  high 
pressure  constitutive  behavior  of  such  a  solid  is  the  plate 
impact  test,  designed  to  induce  this  type  of  planar  shock 
loading.  Both  the  pressure-volume  equation  of  state  (EOS) 
and  the  shear  strength  of  the  solid  can  be  inferred  from 
results  of  this  test,  depending  on  diagnostics  used.  Spall 
fracture  [18]  properties  can  also  be  deduced  if  the 
experiment  is  designed  for  their  interrogation.  The  EOS, 
strength,  and  fracture  properties  all  may  enter  constitutive 
models  for  solids  used  in  hydrocode  simulations  of  terminal 
ballistic  events.  A  comprehensive  and  lengthy  reference 
covering  experiments  and  analysis  of  planar  shock  loading  is 
[19].  The  present  work  places  emphasis  on  brevity  and 
recent  developments  in  theory,  modeling,  and  simulation.  In 
particular,  dynamic  finite  element  simulations  of  shock 
propagation  and/or  spall  fracture  in  polycrystalline 
micro  structures  (e.g.,  [20,  21])  predated  by  [18,  19]  are 
incorporated  herein.  Such  simulations  can  be  used  to  relate 
structure  to  dynamic  properties  important  for  resistance  to 
failure  modes  incurred  during  ballistic  events,  facilitating 
design  of  advanced  material  systems  in  protection  sciences. 

For  each  of  the  two  physical  problems  noted  above, 
classes  of  modeling  technique  are  evaluated  in  terms  of 
flexibility,  complexity,  and  predictive  capability.  Flexibility, 
which  may  also  be  termed  generality,  describes  the  ability  of 
a  model  to  represent  a  breadth  of  physical  behaviors  without 
fundamentally  altering  its  governing  equations.  Complexity, 
i.e.,  sophistication,  of  a  theory  generally  reflects  how 
elaborate  are  the  governing  equations.  Analytical  solutions 
and  numerical  implementation  become  more  challenging  as 
complexity  increases.  Finally,  predictive  capability  refers  to 
50 


a  model's  representation  of  real  physical  behavior  with 
minimal  calibration.  A  particular  theory  is  considered  more 
predictive  than  a  competing  model  if  it  is  able  to  better 
depict  realistic  physics  with  fewer  fitting  parameters  or  ad 
hoc  equations  [22,  23]. 

This  paper  is  organized  as  follows.  Section  2  addresses 
the  ballistic  penetration  mechanics  problem  from  standpoints 
of  classical  hydrodynamic  (1-D)  analysis,  historic  and 
recently  extended  (yet  still  reduced  order,  1  -D)  analysis,  and 
modem  hydrocode  simulations.  Section  3  covers  the  planar 
shock  problem,  including  analytical  solutions  to  the 
governing  equations  for  shock  in  1-D,  steady  wave 
treatments  (also  1-D,  but  usually  requiring  numerical 
solutions),  and  fully  resolved  explicit  dynamic  simulations 
(2-D  or  3-D).  Each  section  logically  progresses  from 
descriptions  of  lower  to  higher  sophistication.  Conclusions 
are  given  in  Section  4,  including  a  table  summarizing  the 
evaluations  of  each  of  the  methods. 

Important  mathematical  relations  pertinent  to  each 
modeling  technique  are  presented  throughout.  Notation  of 
continuum  physics  is  used,  with  vectors  and  tensors  written 
in  bold  italic  font,  and  scalars  and  scalar  components  in  italic 
font.  When  the  index  notation  is  used,  the  normal  summation 
is  implied  for  repeated  indices.  Other  notation  will  be  clear 
from  context. 

2.  BALLISTIC  IMPACT  AND  PENETRATION 

Models  for  the  mechanics  of  terminal  ballistic  events  are 
now  presented.  In  §2.1,  the  classical  1-D  hydrodynamic 
solution  for  steady  penetration  of  a  semi-infinite  body  by  an 
eroding  projectile  is  reviewed.  In  §2.2,  extensions  to  this 
analysis  are  discussed,  most  of  which  remain  1-D  but 
incorporate  finite  strength  of  the  target  and/or  projectile.  The 
modeling  techniques  discussed  in  §  2.1  and  §  2.2  are 
necessarily  restricted  to  the  particular  initial  boundary  value 
problem  of  impact  and  penetration  of  a  rather  long  projectile 
(e.g.,  a  rod  or  shaped  charge  jet)  into  a  deep,  and  therefore 
confined,  target.  In  contrast,  the  computational  modeling 
framework  presented  and  evaluated  in  §2.3  is  capable  of 
describing  impact  events  in  more  complex  physical  systems, 
including  2-D  (often  axisymmefric)  and  3-D  geometries. 

2.1.  Hydrodynamic  Theory:  Classical  Analysis 

The  governing  equation  for  ideal  hydrodynamic 
penetration  of  a  semi-infinite  target  by  a  jet  or  rod  is  now 
derived.  This  derivation  will  be  extended  to  more 
sophisticated  models  in  §2.2. 

The  ideal  hydrodynamic  theory  of  penetration  of  ductile 
targets  by  shaped  charge  jets  was  reported  in  the  early  paper  of 
[24].  The  derivation  rests  on  the  following  primary 
assumptions:  the  penetration  process  is  steady-state  and  one¬ 
dimensional  (1-D);  the  target  is  semi-infinite;  the  projectile  is  a 
continuous  jet;  and  both  target  and  projectile  are 
incompressible  with  null  shear  strength,  i.e.,  are  effectively 
ideal  fluids. 

Denote  spatial  coordinates  at  time  t  by  the  vector-valued 
function  x  =  x(X,  t)  ,  where  reference  coordinates  of  a 
material  point  are  X.  The  particle  velocity  vector  is  the  time 
derivative  of  position: 
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t) 


dx{X,t) 
dt  ' 


(2.1) 

Lo  U~V  Pt 


(2.10) 


Let  a  denote  the  symmetric  Cauchy  stress  tensor,  and 

1 

with  tr(-)  the  trace  operator,  p  =  — -trcr  is  the  Cauchy 

pressure,  positive  in  compression.  The  spatial  mass  density  is 
p(x,t).  In  the  absence  of  body  forces,  the  local  balance  of 
linear  momentum  in  classical  continuum  mechanics  is  [17]. 

V  •  a  =  pi),  (2.2) 

where  the  superposed  dot  is  a  material  time  derivative  and 
where  V(-)  is  the  gradient  with  respect  to  spatial  coordinates. 
The  particle  acceleration  is 

i)(3:,  t)  =  +  Vi)(3:,  t)  ■  i)(3:,  t).  (2.3) 

Under  the  aforementioned  assumption  of  steady  flow, 
v(x,t)  ->  d(:c),  and  the  first  term  on  the  right  side  of  (2.3) 
vanishes,  leading  to 

i)(;c)  =  Vi)(;c)  •  i)(;c).  (2.4) 


The  original  intended  subject  of  this  equation  was 
penetration  of  ductile  metallic  targets  by  ductile  metallic  jets, 
as  in  [24].  Subsequently,  and  often  with  some  degree  of 
success,  (2.9)  and  (2.10)  have  been  invoked  to  describe  the 
steady  penetration  regime  for  long-rod  projectiles  (typically 
of  high  density  and  metallic  origin)  as  well  as  brittle  targets 
(e.g.,  ceramics).  The  advantage  of  this  simple  model  is  that  it 
requires  no  fitting  parameters  whatsoever  for  the  substances: 
only  their  density  ratio  need  be  known.  The  disadvantage  is 
that  it  fails  to  account  for  effects  of  finite  strength  {i.e.,  shear 
stresses)  and  finite  compressibility  of  either  the  projectile  or 
the  target,  both  of  which  may  become  important  depending 
on  the  true  geometry  of  the  target  and  projectile,  the 
materials  involved,  and  the  striking  velocity.  However,  the 
ideal  hydrodynamic  solution  provides  a  useful  limiting  case 
for  comparisons  with  test  data  and  for  comparisons  with 
more  sophisticated  analytical  or  numerical  calculations. 

2.2.  Hydrodynamic  Theory:  Extensions 


For  steady  1-D  flow,  (2.2)  becomes,  with  (t(x)  = 
— (Tj^i(x)  the  axial  stress,  taken  here  as  positive  in 
compression, 

1  dcr  du  j  j  ,, 

- —  =  V  —  ^  —da  =  pvdv.  (2.5) 

p  ax  ax  ' 

Conservation  of  mass  requires 

p  =  pV  -v.  (2.6) 


A  vast  number  of  1-D  analytical  penetration  mechanics 
theories  have  been  based  on  extensions  of  the  ideal 
hydrodynamic  solution  derived  in  §2.1.  Several  notable  such 
models,  termed  reduced  order  theories,  are  now  discussed  in 
this  work. 

Birkhoff  et  al.  [24]  modified  the  hydrodynamic  theory  to 
allow  for  jet  particulation  by  using  a  shape  factor  X  that  has  a 
value  of  unity  for  continuous  jets  and  a  value  of  two  for 
dispersed  particle  jets: 


For  incompressible  flow,  p[:i(X,  t),  t]  =  p[:i(X)]. 

The  physical  problem  is  now  analyzed  by  assuming  that 
an  incompressible  jet  or  rod  with  initial  velocity  V  strikes  the 
target,  an  incompressible  and  infinite  half-space.  The 
stagnation  point  between  projectile  and  target  recedes  with  a 
velocity  of  magnitude  U.  The  axial  stress  P  at  the  stagnation 
point  in  the  projectile  is  obtained  by  integrating  (2.5),  with 
Po  the  constant  projectile  density,  leading  to  [25] 

-  io  do-  =  Po  Iy_u  vdv^  P  =  ipo(F  -  Uf.  (2.7) 


XpoiV  -  uy  =  PtU^  ^  ^  = 
^0 


(2.11) 


Pack  and  Evans  [26,  27]  extended  the  solution  in  (2.10) 
to  allow  for  secondary  penetration,  i.e.,  after- flow  in  addition 
to  primary  penetration,  r.  These  authors  also  included  an 
empirical  correction  for  nonzero  target  strength  Yj : 

l-a^)  +  -.  (2.12) 

PoV^J  Lo  '■  ’ 


Now  considering  the  stagnation  point  in  the  target  with 
mass  density  p-p, 

-  lo  dCT  =  Pt  /y  vdv^  P  =  IptU^-  (2.8) 

Finally,  P  from  (2.7)  is  set  equal  to  P  from  (2.8),  and 
under  the  assumption  of  inviscid  flow  (no  shear  stress),  p  = 
P  .  The  result  is  Bernoulli's  equation  for  steady 
hydrodynamic  penetration: 

p  =  P  =  \po(y-Uf  =\prU\  (2.9) 

The  time  needed  from  first  impact  for  a  projectile  of 
initial  length  Lq  to  fully  erode  is  tg  =  Lo/(U  —  V),  and  the 
final  depth  of  penetration  is  Pq  =  U  •  tg  ■  From  (2.9),  the 
normalized  depth  of  penetration  can  be  obtained  in  terms  of 
the  ratio  of  constant  mass  densities  of  the  target  and 
projectile  materials: 


In  this  theory,  the  scalar  function  a  is  permitted  to 
depend  on  target  and  jet  densities,  and  aYpf(poV^)  =  kR, 
with  k  an  empirical  factor  and  R  the  work  per  unit  volume 
required  for  crater  formation.  Eichelberger  [28]  added  to 
(2.9)  an  empirical  statistical  correction  parameter  y,  and  also 
incorporated  the  net  strength  difference  T/v  =  Yt  ~Yo,  with 
To  the  jet/projectile  strength: 

ypo(V-Uy  =PtU^+2Y^.  (2.13) 

Not  long  thereafter,  Alekseevski  [29]  and  Tate  [30] 
developed  theories  for  long-rod  penetration  of  ductile 
metallic  targets  that  addressed  non-steady  behavior, 
specifically  deceleration  of  the  rod  due  to  finite  strengths  of 
the  projectile  (  Tq  )  and  the  target  ( Rp ).  The  governing 
equation  for  equal  stresses  in  target  and  projectile  at  the 
stagnation  point  is  obtained  in  this  approach  by  modifying 
the  limits  of  integration  in  (2.7)  and  (2.8)  such  that  steady 
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flow  does  not  commence  until  the  stress  reaches  the 
resistance  of  either  material: 

-  ^  dff  =  Po  Cu  ^dv^P=  -  uy  +  Y,;  (2.14) 

do-  =  PtI°  vdv^P  =  ^PtU^  +Rt.  (2.15) 

Axial  stresses  P  are  equated,  which  gives  Tate's  extended 
Bernoulli  equation: 

^Po(V  -uy  +  Yo=^PtU^  +Rt.  (2.16) 


The  complete  theory  developed  in  [29,  30]  includes 
differential  equations  for  projectile  deceleration  and  erosion 
that  must  be  integrated  to  obtain  depth  of  penetration.  Exact 
solutions  have  been  reported  in  [31].  If  deceleration  is 
ignored,  then  the  analytical  solution  for  penetration  depth  is 


Po  _  U  _  1  l'iJ.^lv^+A-V\ 
Lo  U-V  p  \iiV-^V^+a)' 


(2.17) 


where 

P  =  y Pt/Poi  a  =  2(1  —  p.^){Ri-  —  Yq)/Pj.  (2.18) 

The  equivalent  steady  state  solution  can  be  recovered 
from  (2.13)  when  y  ->  1  and  R-p  —  Yg  .  Importantly, 
reduces  to  the  hydrodynamic  result  (2.10)  for  very  high 
velocities  or  low  material  strengths,  i.e.,  when  »  A. 
Analytical  predictions  are  compared  with  numerical  results 
for  metal  rods  penetrating  metal  targets  in  [32].  Such 
comparisons  show  that  target  resistance  Rr  depends  on  the 
experimental  configuration  (e.g.,  geometry)  in  addition  to  the 
target  material's  properties,  implying  that  generally  R-p  ^  Yj-. 

Walker  and  Anderson  [33]  derived  a  time-dependent 
analytical-numerical  technique  to  model  unsteady  long-rod 
penetration  of  semi-infinite  targets.  This  approach  considers 
initial  impact,  requiring  an  initial  interface  velocity  from  the 
shock  jump  conditions  that  will  be  reviewed  later  in  §3.1,  as 
well  as  rod  deceleration.  Key  assumptions  are  invoked 
regarding  the  plastic  flow  field  in  the  target,  obtained  from  a 
dynamic  cavity  expansion  analysis,  and  regarding  the 
velocity  profile  in  the  projectile,  obtained  from  a  priori 
numerical  simulation.  For  a  limiting  case,  the  analogy  of 
Tate's  target  resistance  was  found  to  vary  with  the  dynamic 
ratio  P  of  plastic  zone  size  to  cavity  size: 

Rr=lYr\n[fiy)].  (2.19) 

Dimensional  analysis  of  simulation  results  [34]  for 
metals  showed  that  larger-scale  targets  tend  to  be  weaker 
than  their  small-scale  counterparts  due  to  rate  and  time-to- 
failure  effects.  This  size  effect  arises  since  longer  times  are 
available  for  damage  mechanisms  such  as  cracks  and  shear 
localization  zones  to  incubate  and  propagate  in  larger  targets. 
Similar  size  effects  have  also  been  observed  in  layered 
ceramic-metal  targets  [35].  Experiments  and  hydrocode 
simulations  were  used  to  probe  the  importance  of  strengths 
of  both  target  and  projectile  over  a  range  of  impact  velocities 
in  [36].  Strength  effects  were  found  to  decrease  with 
increasing  impact  velocity,  lending  credibility  to  the  limiting 


Bernoulli  solution  (2.10)  in  the  extreme  hypervelocity 
regime.  Simulation  results  likewise  showed  a  small  effect  of 
compressibility  on  penetration  for  very  ductile  metallic 
targets  [37].  However,  compressibility  becomes  more 
important  in  ceramic  targets  and  concrete,  for  example, 
wherein  initial  porosity  decreases  with  increasingly  high 
compressive  pressure,  and  where  strength  and  pressure  are 
coupled  in  the  constitutive  response  due  to  frictional  effects, 
for  example  [38,  39,  23]. 


Fig.  (1).  Extended  hydrodynamic  theory  based  on  dimensional 
analysis  applied  to  aluminum  oxide  (top)  and  boron  carbide 
(bottom)  [11]. 


Methods  of  dimensional  analysis  have  been  used  to 
provide  further  insight  into  parameters  and  properties 
affecting  the  ballistic  response  of  metallic  [12,  13]  and 
ceramic  [10,  11]  targets.  The  treatment  in  [11]  is  considered 
in  further  detail  in  what  follows  next.  The  penetration  depth 
equation  derived  therein  from  dimensional  and  physical 
considerations  is,  with  aj  fitting  parameters  that  potentially 
depend  on  the  target  material: 


£“  =  i 

t'O  M 

1;  a-^  - 


(ag 


,  V^o/Po  , 


-a]. 


(2.20) 


This  equation  was  applied  to  experimental  penetration 
data  of  [40,  41,  42,  43]  where  target  materials  encompassed 
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the  following  polycrystalline  ceramics:  aluminum  oxide, 
aluminum  nitride,  boron  carbide,  and  silicon  carbide.  All 
experiments  were  conducted  in  reverse  mode  [7],  with  long- 
rod  hypervelocity  impact  and  penetration  into  confined 
cylindrical  ceramic  targets.  Projectiles  were  relatively  pure 
polycrystalline  tungsten.  Impact  velocities  Vg  range  from  1.5 
to  5.0  km/s.  The  primary  discovery  reported  in  [11]  is  that 
normalized  penetration  data  for  all  four  ceramic  target 
materials  can  be  described  well  using  (2.20)  truncated  at 
a-i  =  —a  ,  with  a  single  fitting  parameter,  a  ,  taking  a 
universal  value  of  3.0.  Characteristic  results  are  shown  in 
Fig.  1  for  alumina  and  boron  carbide,  wherein  best  fits  of  a 
provide  little  improvement  over  the  universal  value  of  3. 
With  a  thus  found  independent  of  target  material, 
penetration  depth  depends  only  on  the  properties  of  the 
projectile  (Tq  and  Pn)  and  the  density  of  the  target  entering 
the  factor  Thus,  static  and  dynamic  strength 

properties  and  underlying  mechanisms  would  seem  to  be  of 
little  influence  on  penetration  depth  for  this  test 
configuration,  since  otherwise  a  would  tend  to  vary  among 
the  targeted  materials.  In  future  work,  parametric  numerical 
simulations  of  the  sort  discussed  later  in  §  3.3,  wherein 
material  properties  can  be  varied  systematically  should 
provide  further  insight  into  the  origin  of  a  for  polycrystalline 
ceramics. 

Comparison  of  results  in  [11]  with  another  dimensional 
analysis  [10]  of  a  different  target  configuration  and  velocity 
regime  is  instructive.  This  different  problem  geometry,  as 
explained  in  [7,  44,  45]  for  example,  consists  of  one  or  more 
ceramic  tiles  backed  by  a  semi-infinite  ductile  metal  block. 
Performance  of  the  ceramic  is  measured  by  the  depth  of 
penetration  of  the  projectile  into  the  backing  metal.  In  [10], 
dimensional  analysis  determined  that  penetration  of 
relatively  thin,  metal-backed  ceramic  tiles  could  be  described 
by  two  parameters  that  depend  on  the  type  of  ceramic 
material.  One  parameter  is  needed  to  represent  the  effect  of 
tile  thickness.  Analysis  in  [10]  found  it  to  be  associated  with 
the  ratio  of  fracture  surface  energy  to  elastic  modulus.  The 
second  parameter  describes  the  relationship  between 
penetration  depth  and  impact  velocity  and  appears  to  be 
related  to  the  ratio  of  dynamic  shear  strength  to  target 
density.  Comparison  of  the  dimensional  analyses  in  the  two 
papers  [10,  11]  reveals  an  apparent  transition  from  fracture- 
and  dynamic  strength-controlled  resistance  to  mass  density- 
controlled  resistance  with  increasing  impact  velocity  and 
increasing  target  thickness  or  confinement. 

The  reduced  order  models  discussed  above  in  §  2.2 
present  the  following  positive  features.  Effects  of  strength  of 
the  projectile  and  target  can  be  included  in  the  analysis  with 
relatively  little  increase  in  model  complexity.  Such  effects 
are  often  important  for  lower  striking  velocities,  thinner 
targets,  or  for  stiff  materials  such  as  ceramics.  Solutions  are 
relatively  easy  to  obtain,  either  analytically  or  via  numerical 
quadrature.  Drawbacks  are  that  the  models  remain  one¬ 
dimensional,  such  that  effects  of  lateral  boundaries 
associated  with  finite  sized  projectile  and  target 
configurations  are  omitted.  Furthermore,  parameters  entering 
such  models  are  usually  calibrated  to  match  observed  trends 
or  hydrocode  simulation  results,  rather  than  obtained  from 
first  principles  or  fundamental  experiments  on  constitutive 
behavior.  In  this  sense,  the  reduced  order  models  tend  to  be 


prescriptive  rather  than  predictive,  though  comparison  of 
calibrated  parameters  for  different  materials  subjected  to 
similar  loading  conditions  may  give  qualitative  insight  into 
underlying  physical  phenomena. 

2.3.  Hydrocode  Simulations 

Unlike  the  1-D  models  reported  in  §2.1  and  §2.2, 
computer  simulations  enable  depiction  of  complex 
penetrator-target  geometries,  i.e.,  those  requiring  2-D  or  3-D 
representations.  The  software  used  to  solve  the  governing 
equations  of  continuum  dynamics  is  termed  here  a 
hydrocode,  though  the  materials  involved  need  not  be 
perfectly  hydrodynamic.  In  other  words,  shear  stresses  are 
admitted  in  addition  to  pressure  p,  or  the  stress  tensor  a  need 
not  be  isotropic.  An  early  summary  of  hydrocode 
simulations  can  be  found  in  [9].  The  forthcoming 
presentation  covers  more  advanced  topics  with  a  focus  on 
nonlinear  continuum  mechanics  principles,  including  finite 
deformations  and  nonlinear  material  models. 

The  class  of  models  discussed  next  accounts  for  large 
strains  and  rotations,  as  both  may  occur  during  deformation 
of  ductile  metals  [17,  46],  and  even  in  ceramics  and  minerals 
when  loading  is  predominantly  compressive  [47,  48].  Only 
the  essential  equations  are  provided,  and  Cartesian 
coordinates  are  implied  when  index  notation  is  used.  For  a 
more  comprehensive  treatment  that  encompasses  curvilinear 
coordinates,  see  [17],  with  general  kinematics  addressed  in 
more  detail  in  [49].  Governing  equations  of  finite  anisotropic 
elasticity  are  also  given  in  [50,  51]. 

The  forthcoming  presentation  considers  a 
(poly)crystalline  solid  such  as  a  metal  or  ceramic  material. 
Under  conditions  of  ballistic  impact,  the  material  may 
degrade  in  strength,  i.e.,  undergo  a  damage  process.  In  the 
present  framework,  local  damage  is  represented  by  a  scalar 
state  variable  D  (X,  t)  £  [0, 1] .  Generalization  to  a  vector, 
tensor,  or  multiple  scalars  is  possible  but  not  undertaken  here 
to  maintain  a  terse  review.  The  solid  is  assumed  to  be 
hyperelastic  and  may  undergo  plastic  slip. 

Spatial  coordinates  of  a  deformable  body  are  related  to 
material  coordinates  by  the  time  dependent  motion 

X  =  x(X,t)  =  X +  u(X,t),  (2.21) 

where  u  is  the  displacement  vector.  The  deformation 
gradient  is  decomposed  in  multiplicative  form  as 

F  =  VoX  =  F’^F°F‘’,  (2.22) 

where  F^  includes  thermoelastic  deformation  and 
mechanically  reversible  changes  in  damage  (e.g.,  elastic 
crack  closure  on  load  release),  f  ^  accounts  for  plastic  slip 
from  dislocations,  and  accounts  for  mechanically 
irreversible  damage  mechanisms  like  cracks  and  voids  that 
remain  after  local  elastic  unloading.  A  three  term 
decomposition  of  this  general  form  was  proposed  in  [52]. 
Other  deformation  gradient  representations  that  account  for 
damage  explicitly  include  additive  [53,  54,  55]  and  hybrid 
additive-multiplicative  [56,  57]  forms,  often  derived  or 
motivated  from  homogenization  of  discrete  displacement 
jumps  due  to  subscale  cracks  within  a  volume  element.  The 
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volume  fraction  of  damage  H  is  related  to  the  determinant  of 
F°: 

detF°  =  (1  -  (2.23) 

Besides  its  use  for  solids  with  voids  [58,  52]  or  pores  [38, 
23]  a  multiplicative  damage  term  has  been  introduced  for 
cleavage  cracking  in  crystals  of  metallic,  mineral,  or  ceramic 
origins  [39,  59,  48] 

The  local  balance  laws  of  continuum  mechanics,  in 
spatial  form,  consist  of  the  conservation  of  mass,  linear 
momentum,  and  energy: 

p  =  pV  -v,  V-a+/  =  pv,  pe  =  a-.Vv  —  V  ■  q.  (2.24) 

The  body  force  vector  is  /,  the  internal  energy  per  unit 
mass  is  e,  and  the  spatial  heat  flux  vector  is  q.  Symbols  for 
density,  velocity,  and  stress  are  the  same  as  those  introduced 
in  §2.1.  The  local  balance  of  angular  momentum  is  not 
solved  explicitly,  and  simply  requires  that  the  stress  tensor  a 
be  symmetric.  Point  heat  sources  are  not  relevant  for  the 
present  scope  and  are  omitted  from  the  energy  balance. 
Hydrocode  simulations  integrate  the  coupled  governing 
partial  differential  equations  in  (2.24)  in  time  and  space, 
given  boundary  and  initial  conditions.  The  equations  in 
(2.24)  alone  are  insufficient:  a  constitutive  model  is  required 
for  each  material,  providing  the  mathematical  relationships 
among  stress,  internal  energy,  the  deformation  gradient,  and 
the  latter's  history  and  rate. 

The  thermoelastic  strain  used  in  standard  crystal 
hyperelasticity  [17,  50]  is  the  Green  strain  tensor: 

E  =  i  [(F'^yF'^  -  1],  E,j  =  i  (F^jF^j  -  S,j).  (2.25) 


where  C,jxl...  are  isentropic  elastic  coefficients  of  second- 
and  higher  orders,  Arj  is  entropy  change  measured  from  a 
reference  state  at  temperature  Tg,  and  T/y  are  Griineisen 
coefficients.  Elastic  moduli  depend  on  damage.  The  simplest 
degradation  model  of  the  moduli  is  linear  in  D : 

C[D(X,t),X]  =  [1  -  D(X,t)]Co(X),  (2.29) 

with  Co(X)  =  C(0,X)  the  tensor  of  elastic  moduli  for  the 
undamaged  (poly)crystal  at  the  corresponding  material  point. 
More  sophisticated  approaches  are  needed  to  realistically 
capture  physics  of  arbitrary  loading  cycles,  e.g.,  such  as 
damage  induced  anisotropy  and  differences  in  tensile  and 
compressive  degradation. 

The  viscoplastic  flow  rule  is  often  of  the  general  form 
=  f^(5,'r,{f},f^),  (2.30) 

where  S  can  be  replaced  with  the  most  appropriate  stress 
tensor  for  the  corresponding  configuration  space.  Kinetic 
equations  must  likewise  be  supplied  for  time  rates  of  F^  and 
internal  state  variables,  for  example 

=  F‘’(E,p,{aF‘’l  R)  =  (2.31) 

Dependence  on  elastic  strain,  entropy,  and  internal  state 
variables  is  perhaps  more  physically  meaningfully 
substituted  with  dependence  on  stress,  temperature,  and 
conjugate  thermodynamic  driving  forces.  The  local  balance 
of  energy  in  the  last  equation  in  (2.24),  for  adiabatic 
conditions  often  appropriate  in  impact  dynamics,  can  be 
expressed  as  a  temperature  rate  equation: 

t  =  -rr-.E,  (2.32) 


The  thermoelastic  volume  change  is  measured  by  = 
detf  ^ .  Other  more  recent  formulations  have  used  other 
strain  tensors,  including  Eulerian  and  logarithmic  tensors 
referred  to  material  coordinates,  with  improvements  over  the 
Green  strain  representation  for  metals,  ceramics,  and 
minerals  [47,  60,  48,  61,  62]  Substitution  ofE  with  one  of 
these  alternative  strain  tensors  is  relatively  straightforward. 
Herein,  following  the  standard  approach,  Cauchy  stress  a  is 
related  to  elastic  second  Piola-Kirchhoff  stress  S  via 

a  =  ^F‘^S(E‘^y,  C7ij  =  ^FlkS^.Ff,.  (2.26) 


An  internal  energy  function  per  unit  mass,  with 
corresponding  temperature-entropy  relation,  is 

e  =  e(E,r],{^})  =  e(E,r],pu,D, T  =  de/dr\.  (2.27) 


with  (f }  a  set  of  internal  state  variables  that  affect  the  energy 
stored  in  the  solid,  e.g.,  dislocation  density  and  damage 
D.  This  damage  variable  D  presumably  varies  from  zero  to 
unity  as  the  material  at  the  corresponding  point  loses 
integrity.  The  thermoelastic  stress-strain  relation  for  an 
arbitrary  anisotropic  hyperelastic  solid  is 


^I]  —  P  —  ^IJKL^KL  +  21  ^IJKLMN^KL^MN  + 
'^yiJKLMNPQ^KL^MN^PQ  +  "■  ~  ~  (2.28) 


where  c  is  the  specific  heat  per  unit  volume  at  constant 
thermoelastic  strain,  W  is  the  dissipated  energy  rate  from 
plasticity  and  damage  mechanisms,  and  jS  is  the  fraction  of 
dissipation  converted  to  heat  energy,  i.e.,  the  fraction  of 
stored  energy  of  cold  work  is  1  —  Theory  and  simulations, 
including  a  representative  result  in  Fig.  4  of  [14],  have 
demonstrated  links  among  microstructure  (e.g.,  dislocation 
density),  stored  energy  of  cold  working,  and  ballistic 
penetration  resistance  of  metallic  targets  (see  also  [15]). 

In  addition  to  the  governing  partial  differential  equations 
for  the  bulk  response  of  materials,  physics-based  model 
representations  of  contact  between  impacted  bodies  must  be 
introduced,  accounting  for  momentum  transfer  and  possible 
frictional  effects.  Furthermore,  means  for  addressing 
complete  failure  of  the  material  in  numerical  frameworks 
become  essential  since  element  distortions,  if  too  large, 
render  Lagrangian  finite  elements  inaccurate.  Perhaps  the 
simplest  such  method  is  finite  element  deletion,  where  a 
given  element  is  suitably  eliminated  from  the  calculation 
when  its  strength  vanishes  [63].  This  approach, 
unfortunately,  does  not  always  ensure  that  all  conservation 
laws  of  (2.24)  are  obeyed  consistently  upon  failure.  Another 
approach  involves  conversion  of  Lagrangian  finite  elements 
to  interacting  particles  when  large  strain  thresholds 
associated  with  failure  are  attained  [64].  Results  of  a 
simulation  invoking  this  approach,  specifically  addressing 
dynamic  fragmentation  of  concrete  targets,  are  shown  in  Fig. 
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10  of  [38],  Other  more  recent  methods  such  as  extended 
finite  elements  (X-FEM)  [63,  65]  and  discrete  element 
methods  [66],  which  even  if  less  prevalent  in  hydrocodes, 
may  better  satisfy  the  conservation  laws  and  thus  allow  for  a 
more  realistic  description  of  dynamic  fracture  and 
fragmentation. 

Relative  to  hydrodynamic  theory  and  its  extensions 
discussed  in  §2.1  and  §2.2,  hydrocode  simulations  offer 
many  advantages.  As  noted  already,  complex  geometries  can 
be  resolved,  for  example  heterogeneous,  layered  targets  of 
finite  dimensions  [67,  16].  Stress  wave  interactions  can  be 
visualized  via  modem  post-processing  tools.  Material 
behavior  can  be  sophisticated  and  realistic,  with  nonlinear 
constitutive  models  of  the  sort  discussed  above  enabled.  The 
primary  disadvantage  of  hydrocode  modeling  is  the  cost. 
Firstly,  the  code  or  software  must  be  obtained  or  written,  the 
latter  an  imposing  endeavor  requiring  expertise  in  continuum 
physics  and  computer  science.  Secondly,  solutions  to  most 
impact  problems  are  computationally  expensive,  with  cost 
tending  to  increase  with  model  complexity,  both  geometric 
and  material.  Time  step  size  for  explicit  integration  is 
restricted  by  the  Courant  condition,  which  limits  the  total 
time  duration  over  which  the  momentum  equations  can  be 
integrated.  In  other  words,  impact  events  can  only  be 
simulated  over  short  time  periods.  The  maximum  time  step 
size  becomes  smaller  for  finer  meshes  and  stiffer  materials, 
making  simulations  of  small  and  stiff  structures  very 
expensive.  Finally,  constitutive  models  for  the  material 
response  may  require  calibration,  ideally  via  comparison 
with  data  from  fundamental  experiments  rather  than 
matching  to  ballistic  data.  The  latter  situation  is  sometimes 
unavoidable,  however,  depending  on  availability  of  test  data 
and/or  any  deficiencies  in  the  physical  model  for  the 
material. 

3.  PLANAR  SHOCK  WAVES 

Methods  for  modeling  the  response  of  solids  subjected  to 
shock  wave  propagation,  typically  induced  by  planar  impact, 
are  now  addressed.  In  §3.1,  the  classical  sharp  interface 
treatment  of  a  steady  planar  shock  passing  through  a 
homogeneous  material  is  reviewed,  with  the  Rankine- 
Hugoniot  jump  conditions  the  primary  results.  For  relatively 
simple  kinds  of  constitutive  laws,  the  1  -D  equations  can  be 
solved  simultaneously  and  analytically,  though  not  always  in 
closed  form.  In  §3.2,  another  1-D  treatment  is  presented  for 
analysis  of  planar  shock  waves,  where  steady  state  behavior 
is  presumed,  but  the  shock  width  is  finite  and  is  an  outcome 
of  the  analysis.  The  computational  methods  of  §3.3  resolve 
transient  details  of  the  shock  as  it  travels  through 
heterogeneous  microstructures.  Of  keen  interest  is  the 
predicted  response  of  polycrystalline  materials  pertinent  to 
munitions  (metals)  and  armors  (ceramics)  often  featured  in 
ballistic  applications  of  §2. 

3.1.  Hugoniot  Jump  Conditions  and  Analytical  Solutions 

The  present  analytical  approach  to  modeling  planar 
shocks  in  metals,  ceramics,  and  geologic  materials  involves 
simultaneous  solution  of  the  Rankine -Hugoniot  jump 
conditions  for  conservation  of  mass,  momentum,  and  energy, 
along  with  constitutive  equations  for  the  physical  behavior  of 


the  particular  material  system.  Constitutive  laws  may  be 
elastic,  elastic-plastic,  or  elastic-plastic  with  damage.  The 
onset  of  inelastic  response  corresponds  to  the  Hugoniot 
Elastic  Limit  (HEL);  shocks  of  strength  at  or  exceeding  this 
value  induce  irreversible  deformation  mechanisms  such  as 
slip,  twinning,  pore  collapse,  and/or  fracture. 

For  elastic-plastic  single  crystals,  the  technique  described 
here  and  in  [68]  can  be  applied  only  for  high  symmetry 
orientations:  for  example  shocks  propagating  along  [100] 
and  [111]  directions  in  FCC  and  BCC  crystals  or  along 
[0001]  directions  in  HCP  crystals.  Such  symmetries  reduce 
the  problem  to  simultaneous  solution  of  a  yield  condition 
and  energy  balance  for  the  cumulative  plastic  slip  and 
entropy,  with  the  remaining  conservation  and  constitutive 
laws  sufficient  for  determination  of  the  downstream  material 
state.  Downstream  refers  to  material  behind  the  plastic  shock 
wave,  and  upstream  to  material  ahead  of  the  shock.  For 
lower  symmetry  geometries,  transverse  waves  would  appear, 
and  the  1-D  description  would  be  an  approximation,  often 
severe. 

Considered  is  a  planar  shock  moving  at  natural  velocity 
D  in  the  Lagrangian  direction  X ,  across  which  velocity, 
stress,  and  deformation  gradient  are  discontinuous.  Let  (•)''■ 
and  (•)“  denote  values  of  a  quantity  upstream  and 
downstream  from  the  shock.  The  jump  of  a  quantity  across 
the  shock  plane  is  then 

[[•]]  =  0“  -  0^  (3.1) 

Denote  byv  =  v  —  D  the  velocity  of  the  material  relative 
to  the  shock  front,  with  v  the  particle  velocity  for  1-D 
motion.  The  Cauchy  stress  component  normal  to  the  front  is 
equal  to  the  negative  of  the  shock  pressure:  a  =  —P.  The 
internal  energy  per  unit  mass  is  e,  the  mass  density  p.  Then 
the  Rankine-Hugoniot  equations  for  a  steady  planar  shock 
are  written  as  the  following  jump  conditions  for  mass,  linear 
momentum,  and  energy  [69] : 

[[pvj]  =  0,  [[crj]  -  pv[[v]]  =  0,  [[pv(e  -I-  vV2)  -  crv]]  =  0. 

(3.2) 

The  deformation  gradient  is  uniaxial,  i.e.,  f  =  1  -I- 
(5u/  5A)n  (g)  n,  where  n  is  a  unit  vector  in  the  direction  of 
motion  X  and  u(X,  t)  the  displacement  component  in  this 
direction.  The  above  system  of  equations  becomes  closed 
upon  prescription  of  the  constitutive  model  relating  uniaxial 
deformation  gradient,  internal  energy,  and  stress 
components.  Such  constitutive  models  belong  to  the  general 
class  of  theories  outlined  in  §2.3,  and  for  elastic-plastic 
single  crystals,  those  further  elaborated  in  §3.3. 

Closed  form  analytical  solutions  have  been  derived  for 
ideal  cases.  Perfectly  compressible  fluids  were  analyzed  in 
[70].  A  solution  valid  for  nonlinear  elastic,  anisotropic 
crystals  or  polycrystals  incorporating  the  Green  strain  of 
(2.25)  is  presented  in  [50].  Analytical  solutions  for 
anisotropic  single  crystals  in  the  context  of  (material) 
Eulerian  strain  and  logarithmic  strain  were  first  derived  in 
[47]  and  [48],  respectively.  Subsequent  evaluation  of 
solutions  versus  shock  data  for  a  number  of  single  crystals 
[60,  62]  determined  that  Eulerian  theory  is  often  most 
accurate  for  ductile  metals,  while  logarithmic  theory  is 


55 


Recent  Patents  on  Engineering,  2017,  VoL  11,  No.  1 


J.D.  Clayton 


preferred  for  strong  ceramics  with  a  relatively  high  ratio  of 
shear  to  bulk  modulus. 

Analytical  solutions  can  also  be  obtained,  not  necessarily 
in  closed  form,  for  materials  undergoing  inelastic  behavior 
such  as  plastic  slip.  Isotropic  elastic-plastic  solids  were 
addressed  in  [69,  71].  More  recently,  incremental  analytical 
methods  have  been  developed  to  account  for  multiple  yield 
limits  in  isotropic  polycrystals,  as  shown  in  Fig.  1  of  [72]  for 
the  ceramic  titanium  diboride.  Anisotropic  elastic-plastic 
single  crystals  have  been  addressed  in  [48,  68].  In  this  case, 
the  aforementioned  restrictions  on  symmetry  with  respect  to 
shock  propagation  direction  apply,  and  inelastic  deformation 
is  of  the  form 


equations  relative  to  a  coordinate  frame  that  moves  along 
with  a  steady  {i.e.,  constant  velocity)  wave.  The  steady  wave 
method  has  been  invoked  to  study  plastic  shocks  in  isotropic 
solids  in  [76,  77].  The  first  theory  and  application  of  the 
method  towards  anisotropic  elastic-plastic  crystals  were 
described  in  [78],  with  supplementary  studies  involving  this 
method  as  well  as  finite  difference  simulations  reported  in 
[79,  80,  68]. 


The  system  of  equations  whose  simultaneous  solution  is 
sought  in  a  steady  wave  analysis  consists  of  the  following 
continuity  and  momentum  conservation  laws  [68] : 


du  _  T-v  d/l  dP  _  n 

dr  “  dr '  dr  “  dr' 


(3.4) 


f^(y)  =  exp(y  Xa  <8)  tn“)  =  l  +  r'Za  + 

s‘^  (g)  m“)2  -I-  s“  (g)  m“)^  -b  •••.  (3.3) 

The  cumulative  shear  y  thus  encompasses  the  entire 
deformation  history  in  a  single  kinematic  parameter.  The  slip 
direction  and  slip  plane  normal  vectors  for  system  a  are  s“ 
and  m“.  A  one  parameter  model  for  strength  is  also  typically 
prescribed,  though  history  effects  can  be  admitted  without 
much  more  complexity  [48].  Results  of  this  type  of  analysis 
are  shown  in  Fig.  5  of  [48]  for  single  crystal  diamond 
shocked  along  a  cube  axis,  a  very  strong  solid  which  could 
undergo  slip  or  cleavage  on  its  octahedral  planes  when 
shocked  above  its  HEL. 

Advantages  of  the  method  of  analysis  described  here  in 
§3.1  include  relative  simplicity:  few  material  parameters  are 
needed,  and  solutions  are  obtained  nearly  instantly.  The 
method  is  flexible  enough  to  incorporate  various  nonlinear 
anisotropic  thermoelastic  potentials.  Disadvantages  are  that 
only  isotropic  solids  or  highly  symmetric  orientations  can  be 
modeled  and  time  dependent  effects  such  as  viscosity  (e.g., 
explicit  strain  rate  effects  on  strength)  are  omitted.  Finally, 
because  shock  is  treated  as  a  perfect  jump  discontinuity,  no 
information  regarding  its  structure  such  as  shock  width  or 
values  of  state  variables  within  the  shock  front  is  obtained. 
While  only  a  single  fitting  parameter  may  be  required  as  in 
[68],  its  value  must  still  be  prescribed  via  comparison  with 
shear  strength  data  from  experiments  or  results  of  an 
independent  model  of  material  response. 

For  materials  undergoing  inelastic  response,  this  class  of 
models  can  be  used  effectively  and  efficiently  to  probe 
possible  underlying  mechanisms  such  as  slip,  cleavage, 
twinning,  and  pore  collapse  by  analyzing  various 
orientations  and  comparing  theoretical  yield  criteria  with  test 
data  [73,  74,  75].  The  resulting  information  can  be  used  to 
inform  equations  or  parameters  of  more  sophisticated 
material  models  such  as  those  invoked  in  steady  wave  or 
fully  resolved  finite  element  calculations  of  §3.2  and  §3.3. 
Furthermore,  a  model  can  be  developed  for  one  orientation 
or  symmetry  by  calibrating  to  test  data  then  used  to  predict 
behavior  for  another  configuration  for  which  data  may  be 
unavailable  [23]. 

3.2.  Steady  Wave  Solutions 

The  steady  wave  class  of  modeling  shock  waves  involves 
transformation  of  governing  partial  differential  equations  of 
dynamic  continuum  mechanics  to  ordinary  differential 


The  shock  velocity  is  D,  and  Y  =  X  —  Dt  is  a  coordinate 
moving  with  steady  velocity  in  the  direction  X  of  shock 
propagation.  The  deformation  is  again  presumed  uniaxial, 
i.e.,  F  =  1  +  (du/  dX')n  (g)  n  as  in  §  3.1,  and  1=1-1- 
du/  dX  is  the  axial  compression  ratio.  The  system  of 
equations  also  includes  those  of  the  constitutive  model, 
which  falls  into  the  classes  described  in  §2.3,  and  may  be  of 
any  degree  of  sophistication  so  long  as  symmetry  respects 
the  unaxial  kinematic  condition.  Time  differentiation 
entering  kinetic  equations  for  the  constitutive  model  is 
transformed  to  differentiation  with  respect  to  the  moving 
coordinate  Y.  In  general  the  coupled  system  of  governing 
ordinary  differential  equations  must  be  integrated 
numerically.  This  contrasts  the  less  advanced  material 
models  for  which  analytical  solutions  to  jump  conditions  are 
restricted  in  §3.1.  For  certain  simple  cases  of  isotropic 
material  response,  analytical  solutions  to  the  continuous 
steady  wave  problem  may  be  possible  [76]. 

Advantages  of  the  steady  wave  method  include  the 
following:  a  detailed  description  of  the  steady  shock 
structure  is  obtained,  solutions  are  obtained  at  relatively  low 
computational  cost,  no  artificial  viscosity  is  used  (unlike 
many  finite  difference  and  finite  element  simulations  of 
strong  shocks),  and  sophisticated  rate-  and  temperature- 
dependent  constitutive  models  are  enabled.  An  example  of 
achievable  results  is  shown  in  Fig.  11  of  [78],  which 
demonstrate  effects  of  lattice  orientation  and  shock  pressure 
on  shock  width  in  aluminum.  The  assertion  of  low 
computational  expense  has  been  verified  via  comparison 
with  (more  lengthy)  finite  difference  simulations  of  the  same 
physical  problem  [68].  The  computational  efficiency  of  the 
steady  wave  method  facilitates  its  role  as  a  tool  for 
investigating  and  parameterizing  advanced  constitutive 
models  for  subsequent  use  in  more  expensive,  fully  resolved 
simulation  frameworks. 

Disadvantages  result  from  the  symmetry  and  steadiness 
restrictions.  Accordingly,  effects  of  transverse  waves  for 
non-symmetric  crystal  orientations  are  ignored,  unsteady 
waves  cannot  be  addressed,  and  material  properties  must  be 
spatially  homogeneous.  Only  fully  resolved  simulations  such 
as  finite  difference  calculations  or  explicit  finite  element 
simulations  of  the  class  described  next  in  §3.3  are  capable  of 
quantilying  transient  aspects  of  evolving  shock  waves. 
Important  transient  effects  include  elastic  precursor  decay 
and  wave  interactions,  e.g.,  tensile  reflections  leading  to 
spall  failure.  Fully  resolved  methods  must  be  used  to  account 
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for  heterogeneous  microstructures  and  for  quasi-longitudinal 
and  quasi-transverse  waves  that  arise  in  low  symmetry 
crystal  orientations. 

3.3.  Dynamic  Finite  Element  Analysis  of  Microstructures 

Unlike  the  models  of  planar  shock  propagation  evaluated 
in  §3.1  and  §3.2,  wherein  geometries  of  the  solid  are 
restricted  to  a  single  space  dimension  {i.e.,  1-D),  in  fully 
resolved  dynamic  calculations  of  §3.3,  the  geometry  may  be 
2-D  or  3-D.  Furthermore,  morphology  of  microstructure— 
grain  sizes,  grain  shapes,  grain  boundaries,  secondary  phases 
at  grain  boundaries,  and  so  forth— are  now  resolved 
explicitly.  Geometric  rendering  of  complex  polycrystalline 
geometries  via  finite  elements  is  non-trivial,  with  the 
procedure  described  for  polyhedral  grains,  for  example,  in 
[8 1 ,  82]  and  references  therein. 

The  present  focus  is  on  models  wherein  single  crystals 
within  each  polycrystal  are  resolved  explicitly,  as 
represented  in  crystal  elastic  or  elastic -plastic  simulations  of 
realistic  microstructures  [54,  83,  81,  84,  85,  86].  The 
constitutive  models  for  response  of  the  material  in  a  dynamic 
finite  element  context  are  a  subset  of  the  general  class  of 
nonlinear  continuum  thermomechanics  models  discussed  in 
§3.3.  For  single  crystals  that  undergo  plastic  deformation, 
the  plastic  velocity  gradient  is  the  sum  of  contributions  of 
slip  rates  7“,  where  the  superscript  denotes  a  slip  system  for 
dislocation  glide  with  direction  s“  and  plane  normal  m“: 

iP  =  fPfP-i  =  (8)  m“.  (3.5) 

The  slip  direction  and  slip  plane  normal  are  orthogonal 
and  of  unit  length,  i.e.,  are  those  of  the  crystal  lattice  prior  to 
thermoelastic  deformation.  Denoting  the  resolved  shear 
stress  acting  on  a  system  by  t“,  the  flow  rule  for  slip  rates  is 

Ya  =  Ka(T“,  T,  R}),  t“  =  fa:  [F^5“  (g)  (F^)-Tm“].  (3.6) 

When  regions  of  very  large  defect  density  are 
encountered,  (2.22)  may  be  insufficient  when  is 
attributed  to  dislocation  slip  processes  alone  as  in  (3.5).  An 
intermediate  term,  denoted  here  by  F‘ ,  can  be  used  to 
quantify  residual  lattice  deformation  due  to  defects  within 
the  local  volume  to  which  (2.22)  is  assigned: 

F  =  \IoX  =  F’^F‘fP.  (3.7) 

Here  it  is  assumed  that  damage  is  absent  within  the  single 
crystal;  otherwise,  F®  may  be  appended  to  the  right  side  of 
(3.7).  The  particular  form  of  F‘  depends  on  the  class  of 
defect,  defect  arrangement,  and  scale  of  resolution,  as 
derived  via  a  number  of  theoretical  methods,  [83,  87,  88,  89, 
90,  91,  92,  17].  Analysis  has  demonstrated  the  importance  of 
inclusion  of  F‘  in  the  constitutive  description  when 
dislocation  densities  approach  the  theoretical  maximum, 
which  is  possible  for  severe  plastic  deformation  or  shock 
loading  [73,  91].  Since  is  isochoric  when  attributed  solely 
to  slip,  any  residual  volume  changes  in  the  crystal  are 
omitted  if  not  captured  by  F‘ .  Lattice  rotations  induced  by 
disclinations  may  be  described  by  rotational  part  ofF^  [93], 
and  volume  changes  associated  with  point  defects  such  as 
vacancies  or  interstitial  atoms  may  incorporated  as  well  [94, 
95]. 


Also  not  addressed  explicitly  in  equations  herein  is  the 
possibility  of  deformation  twinning.  This  phenomenon  is 
observed  in  many  kinds  of  crystals  deformed  at  high  loading 
rates,  especially  those  of  lower  symmetry  and  low  stacking 
fault  energy.  Models  invoking  pseudo-slip  kinetics  [96,  97, 
73,  98,  99,  75]  or  phase  field  descriptions  [100,  101,  102, 
103,  104]  have  been  implemented  in  numerical  simulations 
of  twinning  processes  in  single  crystals  and  polycrystals. 
Phase  transformations  can  likewise  be  included  via  diffuse 
interface  modeling  [105]  or  other  prescribed  criteria  in 
continuum  frameworks.  An  example  result  of  a  dynamic 
polycrystal  simulation  incorporating  the  latter  is  shown  in 
Fig.  3  of  [86],  wherein  the  crystal-to-glass  solid-solid 
transformation  in  shocked  boron  carbide  is  modeled  via  an 
intrinsic  nonlinear  elastic  instability  criterion  [106,  107]. 
Quantitative  data  from  3-D  simulations  demonstrating 
validity  of  the  continuum  theory  for  boron  carbide  versus 
planar  impact  experiments  are  shown  in  Fig.  4(a)  of  [86], 
with  further  comparisons  to  various  analytical,  atomic,  and 
experimental  results  listed  in  Table  4  of  that  reference. 

Cohesive  fracture  models  are  now  a  standard 
representation  of  failure  of  the  material  at  the  mesoscale,  i.e., 
the  length  scale  of  polycrystalline  microstructures  of  the 
order  of  the  grain  size.  The  first  explicit  dynamics 
simulations  coupling  finite  crystal  plasticity  with  cohesive 
failure  seem  to  be  those  reported  in  [54,  84]  with  a  follow-up 
study  of  spall  in  [20].  Key  equations,  in  generic  forms,  for 
cohesive  fracture  models  of  the  sort  implemented  in  spall 
fracture  simulations  such  as  that  in  Fig.  5  of  [84]  and  others 
in  [20,  21]  are  reviewed  in  what  follows  next.  Quantitative 
comparison  of  cohesive  finite  element  results  with  spall 
fracture  experiments  is  obtained  by  inspection  of  Figs.  10 
and  26  of  [20]. 

Let  the  crack  opening  displacement  vector  across  two 
crack  faces  initially  coincident  at  point  X  be  defined  as  the 
displacement  jump 

5(A,t)  =  [[u(A,t)]],  (3.8) 

where  displacement  u  need  be  continuous  with  respect  to 
position  only  within  regions  of  the  body  that  have  not 
undergone  fracture.  The  unit  vector  n(X)  is  normal  to  the 
surface  of  impending  fracture.  Let  tg  denote  the  traction 
vector  per  unit  reference  area: 

to  =  P  n,  (3.9) 

with  P  =  (p/pg')aF~^  the  first  Piola-Kirchhoff  stress.  A 
traction-separation  law  is  prescribed  in  the  cohesive  zone  of 
a  generic  functional  form: 

to  =  to(S,{j}).  (3.10) 

Possible  history  effects  are  addressed  by  state  variable(s) 
{x}.  A  magnitude  of  displacement  Sq  is  usually  assigned  as  a 
material  property,  beyond  which  traction  vanishes  and  the 
formerly  cohesive  surfaces  become  free  surfaces.  For 

simulations  invoking  explicit  numerical  integration,  a 

stiffness  matrix  is  usually  not  needed,  in  which  case  the 
traction  function  need  not  have  a  continuous  derivative  with 
respect  to  opening  displacement.  The  work  done  during 
separation  can  be  related  to  a  surface  energy  of  fracture  Y : 
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Y  =  (3.11) 

where  the  path  of  integration  ends  when  the  critical 
separation  magnitude  is  attained.  A  commonly  used  model  is 
the  triangular  cohesive  degradation  function 

ko|  =  c^c(l“|sk^c)  l?ol  =  o'cCl  “  |S|/5c).  (3.12) 

This  law  can  be  invoked  separately  for  each  magnitude  |  •  | 
of  normal  and  shear  components  of  traction  and  opening 
displacement.  The  resolved  scalar  stress  component  at  which 
the  cohesive  zone  starts  to  open  is  Oq,  Le.,  is  the  strength 
required  to  initiate  fracture  or  separation.  Only  two  of  the  three 
parameters  Y ,  Oq  ,  and  Sq  must  be  assigned  since  (3.11) 
eliminates  one  of  these  algebraically.  More  sophisticated 
models  accounting  for  mode  mixity  and  other  physical 
behaviors  have  been  implemented  [108,  109]  but  often  at  the 
cost  of  additional  calibration  or  parameters. 

Cohesive  failure  models  allow  for  realistic  modeling  of 
dynamic  failure  of  microstructures,  including  crack  speeds, 
branching,  and  stress  wave  interactions.  Crack  sizes  and 
shapes  are  fully  resolved.  Positives  of  this  class  of  model 
include  relatively  few  parameters  [minimally  two,  e.g.,  Sq 
and  Oq  in  the  context  of  (3.12)]  and  distinct  behavior  of 
interfaces  and  bulk  material,  meaning  that  traditional  solid 
continuum  elements  can  be  used  for  representing  the  latter. 
However,  more  parameters  are  needed  to  account  for 
distributions  of  strengths  and  surface  energies  among 
potential  failure  sites  in  more  realistic  simulations  of 
heterogeneous  materials.  Anisotropy  and  nonlinearity  of  the 
bulk  crystal  response  can  be  incorporated  via  the  physically 
rigorous  theory  outlined  in  §2.3  and  earlier  in  §3.3. 

Compared  to  analytical  solutions  to  the  Rankine- 
Hugoniot  jump  conditions  of  §3.1,  or  to  numerical  solutions 
to  the  1-D  equations  of  motion  under  the  steady  wave 
assumption  in  §  3.2,  explicit  dynamics  calculations  are 
obtained  at  very  high  expense.  Mesh  generation  procedures 
for  realistic  microstructures  are  often  tedious  and 
cumbersome,  implementation  of  complex  anisotropic  elastic- 
plastic  constitutive  models  is  difficult,  and  integration  of  the 
governing  equations  of  nonlinear  continuum  mechanics  in 
conjunction  with  the  equations  of  these  constitutive  models 
is  computationally  costly.  As  remarked  already  in  §2.3,  the 
Courant  condition  setting  the  maximum  time  step  size 
becomes  prohibitive  for  very  small  element  sizes  which  in 
turn  are  unavoidable  for  modeling  of  polycrystals  of  grain 
dimensions  on  the  order  of  a  few  micrometers  or  less.  On  the 


other  hand,  the  reduced  order  (1-D)  models  in  §3.1  and  §3.2, 
while  enabling  sophisticated  constitutive  models  for  single 
crystals  or  homogenized  polycrystals  with  overall  textures, 
do  not  account  for  details  of  microstructure  morphology. 
Grain  sizes,  grain  shapes,  and  grain  boundary  geometries  are 
all  omitted,  as  are  failure  entities  such  as  discrete  fracture 
planes,  individual  voids,  and  distinct  adiabatic  shear  bands. 
In  contrast,  all  of  these  effects  are  potentially  included  in 
fully  resolved  calculations  of  §3.3.  Thus,  only  the  latter  fully 
resolved  calculations  offer  the  possibility  of  optimization  of 
materials  at  the  micro-  or  mesoscale  for  improved 
performance  (e.g.,  ballistic  penetration  resistance). 
Furthermore,  such  calculations  may  serve  as  a  source  of 
fundamental  functional  forms  and/or  material  parameters  for 
constitutive  model  equations  invoked  at  higher  length  scales. 
See  for  example  [54,  14,  15,  81]  for  examples  of  this 
sequential  multiscale  approach  applied  to  metals  and 
ceramics  for  protection  applications. 

Numerical  implementation  of  cohesive  models  is  deemed 
straightforward,  but  additional  nodal  degrees  of  freedom 
increase  the  computational  expense  where  duplicate  nodes 
are  inserted  along  failure  planes.  A  limitation  is  that  fracture 
paths  are  constrained  to  follow  element  boundaries.  When 
fractures  are  restricted  to  pre-existing  interfaces  such  as 
grain  or  phase  boundaries,  cohesive  finite  elements  can  be 
seeded  a  priori  along  such  pending  failure  surfaces.  Cleavage 
fracture  on  specific  planes  can  also  be  depicted  [110]. 
Extremely  fine  meshes  are  often  necessary  for  representing 
the  small  fracture  process  zones  of  materials  with  high 
strength  and  low  surface  energy;  see  for  example  the  spall 
simulations  of  polycrystalline  silicon  carbide  in  [21]. 
Another  class  of  models  used  to  represent  fracture  of 
microstructures,  not  considered  further  here,  invokes  the 
phase  field  {i.e.,  diffiise  interface)  concept  [111,  112,  113, 
82,  114]  in  contrast  to  the  sharp  interface  assumption 
inherent  in  the  cohesive  finite  element  method. 

4.  CURRENT  AND  FUTURE  DEVELOPMENTS 

Current  modeling  techniques  describing  two  kinds  of 
physical  problem— ballistic  penetration  and  planar  shock 
wave  propagation— have  been  presented  and  compared.  For 
each  kind  of  problem,  models  have  been  categorized  as 
analytical,  reduced  order,  or  fully  resolved.  Perspectives  on 
relative  complexity,  generality,  and  phenomenology  are 
summarized  in  Table  1.  Trends  are  an  increase  in  flexibility 
{i.e.,  increased  fidelity  of  physical  description)  with 
increasing  complexity/sophistication,  albeit  usually  at  the 


Table  1.  Classes  of  impaet  meehanies  models,  examples,  and  general  charaeteristics. 


Analytical  (1-D) 

Reduced  Order  (1-D) 

Fully  Resolved  Dynamics  (2-D  or  3-D) 

Examples:  penetration  mechanics 

Bernoulli  hydrodynamics 

Extended  hydrodynamics 

Hydrocode  simulations 

Examples:  planar  impact 

Hugoniot  solutions 

Steady  wave  analysis 

Explicit  finite  element  simulations 

Complexity  (usual  governing  equations) 

Low 

Moderate 

High 

Generality  (usual  physics  addressed) 

Low 

Moderate 

High 

Phenomenology  (typical  calibration  and  parameters) 

Low 

High 

Moderate 
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expense  of  increasing  number  of  fitted  parameters  or 
equations.  These  trends  are  general,  and  exceptions  of  course 
are  possible. 

It  is  emphasized  that  only  the  fully  resolved  dynamic 
simulations,  though  computationally  costly,  are  able  to 
provide  insight  into  effects  of  fine  scale  microstructures  of 
materials  involved.  Such  effects  are  implicitly  embedded 
within  parameters  entering  reduced  order  and  analytical 
treatments,  but  their  origins  are  not  evident  in  such  cases, 
thus  prohibiting  computationally  driven  materials  design.  In 
contrast,  modem  simulations  at  multiple  length  and  time 
scales  are  becoming  increasingly  important  in  development 
and  optimization  of  advanced  functional  materials  for  high 
rate  applications  such  as  those  related  to  terminal  ballistics. 
Future  developments  will  advance  such  material  system 
design  via  computer  simulations. 
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